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Abstract: The instrument development and design of a prototype frequency- 
domain optical imaging device for breast cancer detection is described in 
detail. This device employs radio-frequency intensity modulated near- 
infrared light to image quantitatively both the scattering and absorption 
coefficients of tissue. The functioning components of the system include a 
laser diode and a photomultiplier tube, which are multiplexed automatically 
through 32 large core fiber optic bundles using high precision linear 
translation stages. Image reconstruction is based on a finite element solution 
of the diffusion equation. This tool for solving the forward problem of 
photon migration is coupled to an iterative optical property estimation 
algorithm, which uses a Levenberg-Marquardt routine with total variation 
minimization. The result of this development is an automated frequency- 
domain optical imager for computed tomography which produces 
quantitatively accurate images of the test phantoms used to date. This paper 
is a description and characterization of an automated frequency-domain 
computed tomography scanner, which is more quantitative than earlier 
systems used in diaphanography because of the combination of intensity 
modulated signal detection and iterative image reconstruction. 
©1997 Optical Society of America 

OCIS code: (170.0170) Medical optics and biotechnology; (170.3010) Image reconstruction 
techniques; (170.3890) Medical optics instrumentation: (170.5280) Photon migration. 
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1. Introduction 


In recent years, there has been a significant amount of research on the development of 
improved methods for optical imaging of human tissues. One of the most promising 
applications of imaging with diffuse light is cancer-detection in the human breast. 
Transillumination, or diaphanography, have been used in various forms for several decades [1- 
4], in an attempt to use the transmitted intensity of red or near-infrared light to image the 
interior of the breast. These early attempts ended with little success mainly because of poor 
image resolution and an inability to detect lesions less than 2 cm in diameter [4-6]. Since the 
design of these early imaging systems, methods to improve image resolution and detectability 
have been realized through the use of computed tomography reconstruction [7-14]. This paper 
describes the essential features of an optical breast cancel' imaging instrument which is being 
developed and will be evaluated in early clinical trials. 

In general, there have been two main theoretical developments in diffuse optical 
imaging, namely improved methods for solving both the forward problem to predict light 
propagation in tissue [15-19], and the inverse problem for estimating the tissue optical 
properties based upon projection measurements along the surface [8-10, 20, 21]. The forward 
prediction of light propagation in a highly scattering medium such as tissue has been modeled 
accurately using several radiation transport methods. Perhaps the most simple and yet 
generally applicable for near-infrared light propagation in tissue is the diffusion model of light 
transport. In the present prototype imaging method, a Unite element solution of the diffusion 
equation has been demonstrated to accurately predict the detected light remitted from tissue 
phantoms, and serves as the basis for the computational reconstruction of images [10-13]. 
Efforts to understand and predict light propagation in tissue, using diffusion theory solutions 
have also been undertaken here. Independent confirmation of how light propagates in tissue 
based on a Monte Carlo statistical approach has been included to confirm the quality of the 
measurements obtained with the new data acquisition system. 

The instrument design has followed many of the recent developments in frequency- 
domain tissue spectroscopy [22-27]. The use of intensity modulated light allows a 
measurement of two parameters, namely intensity and phase, which are influenced uniquely by 
the absorption and diffusion coefficients. These two sets of measurements then permit 
quantitative optical property reconstruction. The major technological objective is to design a 
source-detector array which is low in cost while also being able to acquire multiple 
measurements rapidly. A single source-single detector scheme has been used here which 
incorporates fast multiplexing to minimize the data acquisition time, yet is able to record 256 
projection measurements through the tissue. The system has been specifically designed to 
maximize light detection, thereby maximizing the thickness of tissue which can be imaged. 

This paper outlines the design and development of a fully functional breast cancer 
screening device which should be able to provide quantitative images of absorption and 
diffusion coefficient maps, and the use of multiple projections with a reconstruction algorithm 
should provide increased spatial resolution over diaphanography [10-14]. While this type of 
imaging has been studied extensively both theoretically and in test phantoms, there have been 
few descriptions of automated systems for computational image reconstruction for breast 
cancer detection presented in the published literature [28,29]. This type of system is described 
here along with some representative images obtained from test phantoms. Plans for future 
clinical use are discussed. 

2. Instrument design 

2.1 Hardware 

The light source for the optical imaging system is a diode laser (SDL Inc., CA) at 800 nm 
wavelength and output power of 100 mW. The current to this laser is modulated by mixing 
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continuous and radio frequency (RF) current at a bias Tee (Picosecond Pulse Labs model 
5545, Boulder CO). The RF current is supplied by a signal generator (Marconi Instruments 
model 2023, Ft. Worth TX), typically run at 100 MHz and power output of 13 dBm. Light is 
delivered to and from the tissue phantom via fiber optic bundles with 1/4 inch (approx. 6 mm) 
inner diameter (Cuda Products, FL). 



Fig. 1. Schematic of the automated imaging instrument including hardware and 
software processing. Source optical fibers are indicated in red and detector optical 
fibers in green. 



Fig. 2. Photograph of the optical setup. The 32 fiber bundles set up for taking 
measurements from cylindrical phantoms are connected to two linear translation 
stages one on each side, which are used for rapid multiplexing of light source and 
detector. 
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The light detection is accomplished with a single photomulitiplier tube (PMT) 
(Hamamatsu Inc model R928, Japan), using a specially designed RF modulated housing (ISS 
Instruments, Champaign, IL). The detector is modulated at a slightly offset frequency, so that 
detection can be achieved at the heterodyne beat frequency, which in this case is an offset 
frequency of 100.001 MHz which yields a beat frequency of 1 kHz. 

The source and detector are multiplexed rapidly through the use of 16 source optical 
bundles and 16 detector optical bundles (see Fig. 1 and Fig. 2). The fibers were translated 
laterally in front of both the source and the detector via two linear translation stages, which 
can be moved at 6 inches per second (approx. 15.3 cm/sec) and with a positioning accuracy of 
5 microns (PCLM series SP-4, Anorad Corp.. Hauppage NY). For the detection multiplexing, 
the system has to be able to measure intensity changes over 6 orders of magnitude. To increase 
the dynamic range of intensities that the PMT can detect, an automated filter wheel fitted with 
pre-calibrated neutral density filters was added to the system (Oriel Instruments model 77371, 
Stratford CT). By automatically adjusting the filter before the detector to achieve maximal 
signal, the dynamic range of intensities has been increased from 3 to 8 orders of magnitude. 
An additional advantage of the configuration used here, is that high voltage components are 
spatially separated from the imaging array, which contacts the patient. 

The current signal from the PMT is sampled with a 10 kHz data acquisition board 
using an alternating sampling of the reference and phantom detectors. Before detection, the 
signal is passed through a low pass filter with a cut off of 4 kHz to reduce the high frequency 
noise. 

2.2 Data Aquisition 

All data acquisition, signal processing and instrument control subroutines were 
written within the LABVIEW language (National Instruments, Austin TX). For each 
measurement, a 100 ms sample of both the source and detector signals are taken at 10 kHz, 
and these data samples are processed with software for lock-in detection at 1 kHz to detect the 
real and imaginary components of the AC signal. Preliminary measurements are shown in Fig. 
3 (discussed in Section 2.2) for typical light signals using this short acquisition time. 
Currently the shifting time between samples is 0.5 sec, predominantly due lo the linear 
translation stages. 

Approximately 30 seconds are required additionally for all filter changes during an 
image acquisition. The total data acquisition time for 256 measurements is 7 minutes, at the 
present time. Once all the real and imaginary values are obtained for the 256 source-detector 
pairs, the data is written to a file for image reconstruction. The data aquistion time is longer 
than desired, but should be sufficient to start initial clinical trials for feasability measurements. 
In the future the data aquisition time can be reduced considerably by the addition of more 
photomultiplier tubes for parallel detection. At present all fibers are held in a plastic block 
holder (see Fig. 2), and can be moved manually along the radial direction to adjust to different 
object sizes. The fibers need to be indirect contact with the object being imaged, however 
feasability studies are being done to see if Intralipid filling the space between fiber and tissue 
can be used to help coupling to breast tissue in the clinical tests. 

2.3 Syslem performance 

The performance of the system is dependent upon the characteristics of the modulated light 
source and detector. In particular, a high DC intensity and a high modulation ratio are 
important in order to maintain high signal to noise ratio [27]. The DC intensity can be 
maintained at a high level by simply collecting enough light to achieve a good relative signal, 
however, the modulation ratio is dependent upon the system characteristics. 
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Fig. 3. (a) Measured standard deviation in the AC intensity signal (in volts ) and the 
phase measurement (in degrees) plotted against the total AC intensity (in volts), here 
using a constant DC offset intensity of 1 Volt (5 Volt saturation), and 100 MHz AC 
frequency, and a data acquisition time of 0.2 seconds, (b) Measured modulation 
ratio out of the diode laser as a function of input frequency, using the same 
conditions as in (a). 
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Measurements of AC amplitude and phase error (phase error is in degrees) are plotted 
in Fig. 3(a), as a function of AC amplitude, which demonstrates that in this situation an AC 
intensity of 0.1V and an AC error of 0.001V, results in a S/N ratio of 100, while decreasing 
the AC intensity to 0.01V results in S/N=l. Note here that 0.1 V corresponds to 10% 
modulation of the signal since the DC intensity was fixed at 1.0 V, which is easily achieved in 
imaging through 8 cm tissue phantoms. From these calibration measurements, it is clear that 
maintaining a high AC amplitude at the detector is crucial for reliable data since the S/N ratio 
scales with the square of the AC amplitude. The S/N response is a point of ongoing study to 
better understand the limitations of the system. 

The frequency response of the source-detector system is plotted in Fig. 3(b), 
demonstrating that die optimal frequency range of operation is between 50 and 250 MHz 
where the modulation ratio is approximately 40%. It is interesting to note that if a S/N ratio of 
1 is taken as the useful limit, then reliable measurements can be made with this system up to 
modulation frequencies of 800 MHz, and a S/N ratio of 100 can be maintained up to 500 
MHz. Although it is important to appreciate that at these high frequencies the AC light signal 
is highly attenuated in human tissue, and may not be that useful for imaging. 

Another important parameter to examine is the reproducibility for successive 
measurements of an test object. The two main factors which determine this error are the drift 
of the detectors, and the reproducability of the fiber optic multiplexer. Both of these add to 
the systematic error of the system, however in practice we have found that the drift of the 
photomultiplier tubes is the limiting error presently. The translation stages which are used to 
multiplex the fiber optics have a repeatabilty of approximately 10 microns, allowing very 
small variation in the measured signals. Alternatively the drift in the detectors are difficult to 
quantify exactly but we have found that during the course of repeated measurements on the 
same phantom, that the repeatability is approximately 0.5 degrees in phase and 0.5% in AC 
intensity. This measurement is complicated by the fact that not all measurements have the 
same AC or DC intensities, and that the data aquisition time is flexible. It is important to note 
that these errors are systematic and will manifest themselves in a small shifts of the 
reconstructed optical properties, rather than added noise in the image. 

3. Theoretical tools to analyze photon migration in tissue 

We have developed several numerical tools based on theoretical descriptions of light 
propagation in tissue which have been used to further examine the nature of the data measured 
by our hardware system. In addition, the ability to reconstruct quantitatively accurate images 
depends upon being able to precisely match the measured data with theoretical calculations. 

Theoretically, photon migration in biological tissue can be described by the time 
dependent transport equation [30, 31] 

l ^ r ' f, 'L -V/(r,?,()-ft I(rj,t) + \p(ls')l{r, s< ,t)d£2< + e{r, s,t) > (1) 
c at An j x 

where instantaneous scattering is assumed. The radiance I(r,S, t) is a function of both the 
position in space as well as the direction of the radiated energy within a solid angle dQ. The 
radiance depends on the interaction coefficient u., and the scattering phase function p(s, s' ) . 
The light source is described by £(r , S, t) and c is the velocity of light in the medium. Though 
Eq. (1) provides a valid description for most problems related to optical tomography, there are 
only a few simple cases for which exact solutions can be found. Therefore, the diffusion 
approximation of transport theory is commonly used instead [31]. 


#2827 - $10.00 US 
(C) 1997 OSA 


Received October 8, 1997; Revised December 15, 1997 
22 December 1997 / Vol. 1, No. 13 / OPTICS EXPRESS 397 


While transport theory serves as a reference, the diffusion approximation is more 
attractive from an engineering point of view, because numerical calculations can be performed 
with better efficiency. This becomes particularly important for tomography, where iterative 
methods are used and light propagation has to be simulated multiple times to solve the inverse 
problem. In addition, the diffusion equation can be solved analytically for various geometries 
[32]. Here we briefly summarize their properties, because they complement the experimental 
facilities described in section 2. 



Fig. 4. Cylindrical geometry used for data acquisition: sources and detectors are 
located at the cylindrical boundary within the same plane, A. 

3.1 Computation of radiative transfer by Monte Carlo simulation 

A very flexible way to calculate photon propagation in turbid media is Monte Carlo 
simulation. The probability laws of photon migration can be chosen freely which allows all 
assumptions of transport theory to be incorporated. Photon packages are traced within an 
absorbing and scattering medium and at every interaction the photon weight diminishes 
according to the ratio of absorption and scattering. The scattering characteristic is described by 
the Henyey - Greenstein equation [30], where we typically assume a value of 0.9 for the mean 
cosine of the scattering angle in our simulations. 

The simulation algorithm we use for tracing the photon padi is based on a program 
released previously [33]. The original code was extended to allow the simulation of amplitude 
modulated light sources in the frequency domain. To this end we incorporated a similar 
approach to the one outlined in reference [34]: the propagation distance along the photon path 
is accumulated. At the detector site the path lengths can be associated with a phase delay. 
Phase and photon weight define a complex value. The interference of all contributing photon 
packages then provides an estimate of the detected intensity. 
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The simulation program also allows the definition of different kinds of geometries, 
including a cylinder with absorbing boundaries (see Fig. 4) which corresponds to the phantom 
we use in our imaging experiments. 

3.2 Solving the diffusion equation 

The diffusion approximation turns the problem of photon migration into a simple boundary 
value problem. For a harmonically modulated light source we seek a solution of [31] 

V ■ [flV^r, a>)\ - \p a + — j 0, co) = - S(f, co) (2) 

with boundary conditions of type 

<t> = oc^ (3) 


ac amplitude 



0 20 40 60 80 100 120 140 160 180 



0 20 40 60 80 100 120 140 160 180 
angle source - detector [dee] 

(b) 

Fig. 5. Comparison between experimental data, Monte Carlo simulation and finite 
element calculation for a homogeneous cylindrical phantom (absorption 
u, a ^) .003mm" 1 , diffusion coefficient D=0.72mm). AC amplitude (a) and phase (b) of 
the measurement correspond fairly well to both theoretical models. 
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where n is the normal of the boundary surface. The solution is given in terms of the fluence 
())(?, CO ) which depends on the modulation frequency to, the absorption coefficient ^ a , and the 
diffusion coefficient D. The light source is given by S(r,Co) . The coefficient a describes the 
structure of the boundary. 

To solve Eq. (2) numerically we use a finite element (FEM) algorithm. This 
algorithm calculates the forward problem of photon migration, within the diffusion 
approximation. It is also incorporated into the iterative reconstruction method which retrieves 
the inhomogeneous distributions of the absorption and scattering coefficient by minimizing the 
difference between simulated and measured intensity at the detector locations as well as the 
total-variation over the entire image. This minimization is performed by a Levenberg- 
Marquardt algorithm. Details of the finite element algorithm as well as of the reconstruction 
method can be found in references [10-13]. 

Fig. 5 shows calculations for a homogenous cylindrical phantom. On the one hand, 
these calculation prove, that our theoretical assumptions, incorporated into our Monte Carlo 
simulation meet the experimental conditions. On the other hand, they demonstrate, that the 
diffusion approximation provided by the FEM calculation is indeed a very good 
approximation for the problem with which we are concerned. This is even true at the 
boundary, where the detectors are located, although in general the diffusion approximation is 
only valid far from sources and boundaries [31]. The measured data here fits the FEM 
calculation better than the Monte Carlo simulation, however the differences in both cases are 
small and is likely attributable to systematic problems with matching the boundary conditions 
of our tissue-simulating phantom. 

4. Image reconstruction 

The FEM routine for image reconstruction is called from within LABVIEW and 
processes the data into images of absorption and scattering coefficients. The reconstruction 
program works to minimize both the least-squares difference between the measured data and 
the calculated data from a 2-dimensional FEM of the diffusion equation, as well as minimizing 
the total variation of the image [10-13]. In practice the data fitted is the real and imaginary 
components of the AC signal, which is equivalent to fitting to the AC magnitude and phase. In 
these measurements, there was an inter-fiber intensity variation which could corrupt the data. 
For this reason, all of the following reconstructions have used data measurements which have 
been normalized to the homogenous sample. In practice, normalized data is created by 
precalibration measurements on a homogenous phantom which closely mimics breast tissue. 
In the future, this may also be achieved by imaging the same object at two different 
wavelengths [29], and normalizing the measurements. The images shown below were taken 
with an early version of the system which employed manual multiplexing of the source and 
detector. The computational time per iteration was approximately 2 minutes, but this can be 
decreased to approximately 1 minute per iteration through the use of a dual mesh 
reconstruction algorithm [10]. A total of 10 iterations were obtained for each of the images 
reported here. A minimization method described previously [13] was used with total-variation 
and regularization parameters of 10" 7 and 10, respectively. 

Fig. 6 shows the reconstruction from a cylindrical 25 mm diameter high contrast 
absorbing inhomogeneity placed in the center of a 86 mm diameter scattering phantom. The 
background optical properties were D = 0.67 mm and u. a = 0.0062 mm" 1 , with an 
inhomogeneity having D = 0.67 mm and u. a = 0.056 mm" 1 . In this case the reconstruction is 
quantitatively correct for both absorption and diffusion coefficient, and the full width at half 
maximum for the inhomogeneity is 22 mm. Fig. 7 shows the reconstruction from the same 
size object as in Fig. 6, but with optical properties of D=0.2 mm and u. a = 0.0062 mm" 1 . In 
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Fig. 7. Reconstruction of a 25 mm diameter scattering object (corresponding to a 
decrease in diffusion), with the same background optical properties as in Fig. 6 but 
with inhomogeneity properties of D = 0.2 mm" 1 , and m = 0.0062 mm" 1 . 

5. Discussion 

In the past 7 or 8 years, there have been many reports published discussing the possibilities of 
near infrared tissue imaging, yet very few papers have addressed the technical challenges of 
developing an integrated system for this purpose [28, 29]. One of the obvious applications for 
near infrared computed tomography is breast cancer detection, mainly because recent 
developments in diffuse computed tomography may afford better means of tissue 
discrimination over older clinical devices. Unfortunately, better methods for image 
reconstruction can also be more sensitive to measurement error, and it is crucial to design an 
optical system in which both the random errors of the measurements and the systematic 
mismatch between measurement and theory are all minimized. In this paper, these errors are 
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characterized for our prototype system, and preliminary images of tissue phantoms have been 
presented. 

Systematic errors in optical signal detection need to be eliminated by careful 
normalization of the data or by some other means (i.e. multiple wavelengths, frequencies, etc.). 
In this multiple source-detector arrangement variations in transmitted intensity between fibers 
ranged up to 30%. This level of systematic intensity variation must be corrected in the 
measurements before the data can be used in the reconstruction algorithm. In the cases 
presented here the heterogeneous phantom image data was normalized by data collected from 
a homogenous phantom. At this point, it is still not clear whether the limiting noise in the 
system will be systematic or statistical, although this depends upon the signal from each 
source-detector pair. At detector locations near the source, systematic error tends to impose 
the limitation on precision, while at detectors far from the source, where there is not much 
light signal, statistical errors from the detector dominate. In spite of these limitations, accurate 
absorption coefficient images can be created with the system for imaging phantoms which are 
similar to breast tissue [36-38]. 

The reconstruction of the absorbing object is quantitatively correct and is 
encouraging for images based upon changes in absorption (see Fig. 7). The second image of a 
diffusing object does not reconstruct properly, because there is a slight artifactual 
reconstruction of the diffusing effect into the absorption coefficient image (i.e. the presence of 
a pure absorbing and diffusing inhomogeneity cannot be exactly separated given the 
measurements obtained with this instrument, see Fig. 7). These results suggest that 
measurements are required which better discriminate between absorption and scattering, such 
as the use of higher frequency modulation [35]. Another possibility might be to increase the 
total number of independent measurements through the use of more source-detector pairs. 
Alternatively, biasing the regularization method to better fit the phase shift can preferentially 
enhance the contrast of a diffusing object over an absorbing object [35], and this type of 
approach may turn out to be the most practical. While the inability to quantitatively 
reconstruct diffusion coefficient changes is problematic, most changes between normal and 
tumor tissue are thought to be predominately from blood volume [26, 36], and hence changes 
in absorption would dominate in this case. It is important to note that while the images shown 
in this paper were taken with manual multiplexing of the source and detector fibers, that 
identical images can be obtained presently with the automated system. Also there has not been 
a conclusive test as to how well the reconstruction algorithm will work in highly heterogenous 
media such as real breast tissue, however this will be one of the major goals of the early 
clinical testing. 

In summary, the system presented here is a complete description of a computed 
tomography system which will be used in clinical trials to evaluate the potential for breast 
cancer screening. The architecture of this system provides an automated and accurate method 
for data collection and image generation. 
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